Quantum Monte Carlo calculations of symmetric nuclear matter 
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We present an accurate numerical study of the equation of state of nuclear matter based on 
realistic nucleon-nucleon interactions by means of Auxiliary Field Diffusion Monte Carlo (AFDMC) 
calculations. The AFDMC method samples the spin and isospin degrees of freedom allowing for 
quantum simulations of large nucleonic systems and represents an important step forward towards 
a quantitative understanding of problems in nuclear structure and astrophysics. 
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The equation of state (EOS) of nuclear matter repre- 
sents a challenge in both nuclear structure physics and 
astrophysics. The knowledge of the properties of nuclear 
matter, and in particular of asymmetric nuclear matter, 
is needed to predict the structure, the dynamics and the 
evolution of stars, in particular during their last stages, 
when they become ultra-dense neutron stars. Depending 
on the EOS, the density of nuclear matter in the inner 
shells can reach up to 9 times the core density of stable 
nuclei, po = 0.16 fm^'^ 1]. 

One important step towards the understanding of these 
astrophysical problems is the study of the symmetric nu- 
clear matter, related to the various model NN interac- 
tions available. While considerable advances have been 
made;2, ,3|], it is still impossible to firmly ascertain the 
degree of accuracy of the approximations one has to in- 
troduce in the many-body theories, and substantial dis- 
crepancies still exist among the different theoretical es- 
timates of the EOS, the response functions and Green's 
Functions of nuclear matter. 

Experimental data on symmetric nuclear matter are 
limited to the volume and the symmetry energy of the 
Weizsacker mass formula, and to the nuclear matter 
compressibility. Instead, an indirect test for the the- 
oretical predictions of the EOS of asymmetric nuclear 
matter is provided by the mass-radius relation of a 
neutron star[^ Q, obtained by solving the Tolman- 
Oppenheimer-Volkov equation. 

At present the theoretical uncertainties on the equa- 
tion of state, coming from the approximations one has 
to introduce in the many-body methods, and the lack of 
knowledge of the nuclear interaction, do not allow for def- 
inite conclusions when comparing with astronomical ob- 
servations. However, the recent success in predicting the 
properties of light nuclei gives us some confidence that 
the non-relativistic description of nuclear matter based 
on effective potentials fitted to reproduce NN data and 
the binding energy of light nuclei can be reliable enough. 
The main feature of such nucleon-nucleon interactions, 
besides the short range repulsion, is the explicit depen- 
dence on the relative quantum state of the nucleons which 
can be described using spin and isospin, angular momen- 



tum, spin-orbit, and tensor operatorsQ. 

Properties of light nuclei (A<6) can be efficiently 
computed with high accuracy using modern few-body 
techniques 0, H, 0] or with the ab initio no-core nuclear 
shell model (with A<12[lol|). Quantum Monte Carlo 
techniques based on recasting the Schroedinger equa- 
tion into a diffusion equation (Diffusion Monte Carlo or 
Green's Function Monte Carlo), allowed for performing 
calculations up to A<12[lll [l^]- However, the compu- 
tational resources needed for such simulations are very 
large, because of the summation over all the possible 
states necessary to evaluate the terms of the Hamilto- 
nian with a quadratic dependence on spin and isospin. 
The number of such terms, and the CPU time needed to 
calculate the m, g rows exponentially with the number of 
nucleons; ^^c[lr| or 14 neutrons is the limit for the 
currently available computational resources. 

Since the spatial degrees of freedom are already sam- 
pled, one would like to replace the sum over the spin- 
isospin states with an efficient sampling method. The 
simulation of symmetric nuclear matter requires a mini- 
mum of 28 nucleons in a box replicated in space (7 nucle- 
ons for each spin-isospin state) to obtain a wave function 
with closed shells of momenta, and this is out of the reach 
of the standard Quantum Monte Carlo methods. 

In this paper we show that the spin-isospin is efficiently 
sampled by using the Auxiliary Field Diffusion Monte 
Carlo method ij] , which is based on the use of auxiliary 
variables to linearize the quadratic spin-isospin operators 
of the nuclear matter Hamiltonian, making them treat- 
able in a diffusion Monte Carlo scheme. Up to now it has 
been applied to simulate pure neutron matter (up to 114 
neutrons) jlSt 116], and neutron dropsjlTl Il8l| interacting 
with realistic two- plus three-body interactions. 

Here we extend calculations to include isospin degrees 
of freedom and to deal with the strong tensor-isospin 
force, responsible of the nuclear binding. The method can 
readily handle an asymmetry in the number of neutron 
and protons or the deformation of heavy nuclei. 

In this letter we show that simulations of symmetri- 
cal nuclear matter interacting via a semiphenomenolog- 
ical two-body interaction including spin-isospin depen- 
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dent and tensor components have led to an EOS which 
shows significant difl^erences with respect to that obtained 
within Fermi Hypernetted Chain and Brueckner Hartree 
Fock methods[19], particularly at high densities. Even 
more important is the finding that Quantum Monte Carlo 
simulations do not lead to any lowering of the FHNC or 
BHF energies at p ^ pq. This fact points toward an 
inadequacy of commonly used three-nucleon interaction 
models in the whole range of density. 



Auxihary Field Diffusion Monte Carlo (AFDMC)[l4| 
is an extension of the standard Diffusion Monte Carlo 
method in which the ground state of an Hamiltonian 
H is obtained by solving the imaginary time dependent 
Schroedinger equation 



(1) 



The solution is obtained by evolving a population of 
configurations of the system ("walkers") X = {i?, a, r}, 
where R — {ri , . . . , r/v , } , a ~ {(Ji , . . . , (Tjv } , and t = 
{fi, . . . ,fAr}, with F{X,t) = «'T(^)*(X,t), according 
to 



F{X,t) 



dX'§^^Go{X,X',t)F{X',0) (2) 



The function is a "trial" wave function, usually deter- 
mined by means of variational calculations, and Go is an 
approximation to the Green's function of the imaginary 
time Schroedinger equation: 

(3) 

where D = /2m, Eq is an estimate of the ground 
state energy of the system, and V{X) is the nucleon- 
nucleon interaction. For a long enough imaginary time 
the distribution of the walkers converges to the prod- 
uct ^TiX)'ifo{X) where ^'o is the wave function of the 
ground state of H. This fact allows the computation 
of matrix elements (\I't|0|^'o) of observables O of in- 
terest in a Monte Carlo way. When O = H the value 
obtained is the exact ground state energy of the system. 
The presence of an interaction V(X) including operators 



like (3(Ti • TijBj 



Si ■ Sj) and • Tj is the origin of 



the computational cost in the standard approaches. The 
spin-isospin dependent part oiV{X) {Vsid) can be writ- 
ten as a sum of a matrix ^m.^/j multiplied by spin-isospin 
operators as follow: 



Vsid 



1 



1 



2 ^» 



2 ^ E ^-A„, (4) 

a— 1 n— 1 



where An are the eigenvalues obtained by diagonalizing 
the matrix A, and Sna are operators written in terms of 
eigenvectors of A as follow: 



(5) 



AFDMC uses the Hubbard-Stratonovich method to 
transform the operators S which are quadratic in the 
spin and isospin into linear operators: 



-(l/2)iAS^ _ 



1 



dye 



-\ts 
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Then S are operators which are linear combinations of 
the spin and isospin operators for each nucleon, and A de- 
pend on the interaction. The transformed Green's Func- 
tion is applied to the spin-isospin part of the wave func- 
tion, and its effect consists of a rotation of the spin and 
isospin degrees of freedom (written as four-component 
spinors in the proton-neutron up-down basis) by a quan- 
tity that depends on the auxiliary variable y along with 
multiplication of the state by an overall factor. The sum 
over spin and isospin is replaced by sampling a set of 
rotations of the variables. This procedure reduces the 
dependence of the computational time on the number of 
nucleons necessary for performing a simulation step from 
exponential to cubic. It is therefore possible to perform 
on a regular workstation or on a modest PC cluster cal- 
culations that would require Tflop supercomputers with 
the standard methods. This method, like other diffusion 
Monte Carlo methods, suffers from the so-called "sign 
problem" when it is applied to fermions, and when com- 
plex wave functions need to be used. In our calculations 
we apply the fixed-phase approximation to overcome this 
problem [i^. While this method has already been suc- 
cessfully applied to pure neutron matter [lij. it has not 
been previously used for mixed proton and neutron sys- 
tems. It should be noted that it does not guarantee an 
upperbound to the mixed energy used here. As a test 
for the correctness and the efficiency of our approach we 
reproduced within 0.3 MeV total energy the existing re- 
sults for the binding energy of *He with potentials of the 
vq class [21|. We have also been able to compute binding 



energies for ^^O with this method. 

A crucial point in dealing with nuclear matter is 
the choice of the interaction among nucleons. As al- 
ready mentioned, several modern two-body potentials 
are available nowadays, all fitting the NN data with 
~ 1. We use the potentials of the Argonne class with 
n operators ( AVn) [22l| . While the full version contains 
n = 18 operators, most of the physics is reasonably well 
described by the first 6 operators made up of 4 central 
spin-isospin dependent components and two tensor ones, 
which include the long range one-pion-exchange force. 
The most important missing terms are the spin-orbit 
components. In nuclei and neutron drops the spin-orbit 
contribution to the energy amounts to a few tenths of 
MeV/nucleon. A correction of the order IMeV/nucleon 
can be attributed at low densities to the remaining terms 
included in AV18. Specifically, we have used the interac- 
tion Argonne v'g [2l| truncated by dropping the spin-orbit 
terms, and including only the first six operators, which 
we denote as "our AV6'" . 



TABLE I: AFDMC energies per particle in MeV of 28, 76 and 
108 nucleons in a periodic box at various densities. 





E/A(28) 


E/A(76) 


E/A(108) 


-10 


0.5 


-7.64(3) 


-7.7(1) 


-7.45(2) 


> 


3.0 


-10.6(1) 


-10.7(6) 


-10.8(1) 
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It is well known that two-body NN interaction under- 
bind light nuclei, and one needs to add a specific efi^ective 
three-body potential to reproduce their low energy prop- 
erties. Semi-phenomenological three-nucleon interac- 
tions following the lowest order three-nucleon diagrams 
with one and two intermediate Delta resonance states 
provide a very satisfactory description of the ground state 
energy and the low level spectra of light nuclei up to 
^^C[ll|. We have disregarded such three-body forces in 
our simulations. In nuclear matter they are essential to 
reproduce the experimental saturation density, and, in 
general, they contribute about 10% of the total binding 
energy. A full comparison with the available experimen- 
tal data goes beyond the scopes of the present paper. 
Here we are interested in showing the efficiency of the 
AFDMC methods in dealing with nuclear matter models 
which include realistic tensor interactions like in our AV6' 
potential. Nuclear matter calculations with Argonne Wg 
and Urbana three-nucleon interaction are in progress. 

The results of the calculations with A=28 include box 
corrections that have computed by adding to the two 
body sums contribution of nucleons in the first shell of 
periodic cells. Such procedure is effective. In order to 
assess the magnitude of finite size effects we performed 
calculations with 76 and 108 nucleons at densities p = 
0.08 fm~^ and p — 0.48 fm~^. Results are shown in tabic 
m As it can be seen the results coincide with the ones 
obtained with 28 nucleons within 3 percent. 

In the case of 28 nucleons for each density we gen- 
erated and then propagated a set of 1000 walkers for 
different time-steps ranging from At = 5 x 10~^MeV^^ 
to At = 2.5 X 10~^MeV~^. Each propagation at each 
time-step were performed up to at least a total imagi- 
nary time of t = 2MeV"^ The AFDMC energy is de- 
termined by extrapolating to At —^ 0. In order to lower 
statistical errors, in some case longer total propagation 
time was needed, up to a maximum of i = 6MeV~^ in 
particular at higher densities. Using a parallel supercom- 
puter (typically 16 CPU are employed) a propagation of 
20000 steps requires about 80 processor hours. Then for 
a fixed density we estimated that a maximum of 5000 
CPU hours are needed. In the case of 76 and 108 nu- 
cleons we performed calculations only at a one time-step 
At = 10~^MeV~^ and we propagated until a total imag- 
inary time oi t — lMeV~^. 

We computed the EOS of symmetric nuclear matter 
in the range of densities 0.5 < (p/po) < 3, and com- 




FIG. 1: (color online). Equation of state of symmetric nu- 
clear matter calculated with different methods. Red cir- 
cles represent AFDMC results with statistical error bars and 
the green line is the fitted functional form described in the 
text. Dashed lines correspond to calculations performed with 
other methods ^19!] (blue line with squares; FHNC/SOC; ma- 
genta with diamonds: BHF). Blue triangles represent the 
FHNC/SOC energies corrected by including the low order 
of elementary diagrams as described in the text. Blue ar- 
rows show the corresponding energy shift, which increases at 
higher densities. 



pared it with previous available results obtained with 
the same potential using Fermi Hypernetted Chain in 
the Single Operator Chain approximation (FHNC/SOC) 
and the Brueckner-Hartree-Fock (BHF) in the two-hole 
line approximation To']. AFDMC calculations were per- 
formed with 28 nucleons, filling the shell of plane waves 
with momentum of modulus 1 and providing a wave func- 
tion yielding an isotropic density. 

The results are summarized in Fig. 1 and reported 
in Table [Hi The comparison of the various EOS sug- 
gests the following comments: FHNC/SOC leads to an 
overbinding at high density. A similar indication was 
found by Moroni et al.[23[ after a DMC calculation of 
the EOS of normal liquid "^He at zero temperature, with 
a guiding function including triplet and backflow correla- 
tions. The comparison with the equivalent FHNC/SOC 
calculations of refs 24 , [2F 
ancies 



have shown similar discrep- 
There are two main intrinsic approximations in 
variational FHNC/SOC calculations, which violate the 
variational principle. The first one consists in neglecting 
a whole class of cluster diagrams, the so called elemen- 
tary diagrams, which cannot be summed up by means of 
FHNC integral equations. We have calculated the low- 
est order diagram of this class, namely the one having 
only one correlation bond and four exchange bonds. The 
results obtained show a substantial effect from this dia- 
gram and bring the FHNC/SOC estimates very close to 
AFDMC results, as shown in Fig. 1. The second approx- 
imation is related to the non-commutativity of the corre- 
lation operators entering the variational wave function. 
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TABLE II; AFDMC energies per particle in MeV of 28 nucle- 
ons in a periodic box at various densities. 



p/po 


0.5 


0.75 


1.0 


1.25 


1.5 


1.75 


E/A 


-7.64(3) 


-9.81(4) 


-11.5(1) 


-13.0(1) 


-13.73(7) 


-14.1(2) 


p/po 


2.0 


2.25 


2.5 


2.75 


3.0 




E/A 


-14.0(3) 


-13.5(3) 


-12.7(2) 


-11.7(2) 


-10.6(1) 
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The only class of cluster diagrams contributing to such 
non-commuting terms, which can be realistically calcu- 
lated, is that characterized by single operator chains. It 
is believed that such an approximation is reliable in nu- 
clear matter, but there is no clear proof of this. 

BHF calculations of ref.[l3| predict an EOS with a 
shallower binding than the AFDMC one. It has been 
shown for symmetric nuclear matter, using the AV18 and 
AV14 potentials, that contributions from three hole-line 
diagrams add a repulsive contribution up to ~ 3MeV 
at densities below /h i f2(]| . and decrease the energy at 
high densities [2^. Such corrections, if computed with 
our AV6' potential, would probably preserve the same 
general behavior, and bring the BHF EOS closer to the 
AFDMC one. Therefore, our calculations show that the 
two hole-line approximation used in Ref. is too poor, 
particularly at high density. 

The AFDMC equation of state was fitted with the fol- 
lowing functional form: 

^ = ^+a{x~xf+(3{x^x)\ (7) 

where x — p/ po and the various coefficients are given by 
Eq/A = -14.04(4) MeV, a = 3.09(6) MeV, (3 = -0.44(8) 
MeV, and x = 1.83(1). The resulting compressibility 
K = {E/A) /dx^)- at saturation density x is ~ 

190 MeV. The fit of the EOS allows for computing the 
pressure vs. density for symmetric nuclear matter. 

The availability of an efficient and relatively fast pro- 
jection algorithm for the computation of energies and 
other observables of dense hadronic matter enables the 
possibility of a more quantitative understanding of the 
properties of neutron stars and supernovae, as well as 
that of medium-heavy nuclei. Computations on such sys- 
tems are at present out of reach of the standard GFMC 
methods and available supercomputers. Therefore, the 
extension of AFDMC algorithm to deal with nuclear mat- 
ter is a significant step forward. Some technical im- 
provements on the calculations presented here, such as 
the addition to our AV6' of spin-orbit terms and three- 
body interactions are already underway. The treatment 
of asymmetric nuclear matter, particularly important for 
the determination of the properties of neutron stars, is 
also straightforward, and will be the subject of future 
exploration. 
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